clear; clc; close all;
% 绘制 PTM 随参数 c 变化的分岔图和 Lyapunov 指数

% PTM参数
a = 0.3;
b = 1;
c_values = -4:0.005:1;  % c的变化范围
x0 = [0, 2];  % 初始值数组
N = 1000;  % 计算Lyapunov指数的步数
transient = 200;  % 舍弃的暂态步数
iterations = 1000;  % 用于绘制分岔图的迭代次数

% 存储Lyapunov指数
LEs = zeros(length(x0), length(c_values));  % 初始化Lyapunov指数矩阵
% 存储分岔图
bifurcation_data_0 = zeros(length(c_values), iterations);  % 对于x0 = 0的分岔图
bifurcation_data_2 = zeros(length(c_values), iterations);  % 对于x0 = 2的分岔图

% 对每个初始值进行计算
for i = 1:length(x0)
    for j = 1:length(c_values)
        c = c_values(j);
        
        % 计算Lyapunov指数
        LEs(i, j) = LEs_ptm(x0(i), a, b, c, N, transient);  % 根据初始值x0(i)计算Lyapunov指数
        
        % 绘制分岔图
        x = x0(i);  % 初始值
        for k = 1:transient  % 舍弃暂态
            x = a * (x + b * sin(3 * x) + c)^2;
        end
        
        % 记录后续的值用于绘制分岔图
        for k = 1:iterations
            x = a * (x + b * sin(3 * x) + c)^2;  % PTM更新
            if i == 1
                bifurcation_data_0(j, k) = x;  % 保存x0=0的分岔点
            else
                bifurcation_data_2(j, k) = x;  % 保存x0=2的分岔点
            end
        end
    end
end

% 绘制Lyapunov指数图
figure;
subplot(2, 1, 1);  % 上图为Lyapunov指数图
plot(c_values, LEs(1, :), 'r', 'LineWidth', 2); % 初始值为0的Lyapunov指数
hold on;
plot(c_values, LEs(2, :), 'b', 'LineWidth', 2); % 初始值为2的Lyapunov指数
yline(0, 'k--', 'LineWidth', 1);  % 零刻度线
xlabel('c');
ylabel('LEs');
title('李指数随c变化');
legend({'x_0 = 0', 'x_0 = 2'}, 'Location', 'best');
grid on;

% 绘制分岔图
subplot(2, 1, 2);  % 下图为分岔图
hold on;

% 绘制x0 = 0的分岔图（红色）
for j = 1:length(c_values)
    plot(c_values(j) * ones(1, iterations), bifurcation_data_0(j, :), '.', 'Color', 'r', 'MarkerSize', 1); % 红色分岔点
end

% 绘制x0 = 2的分岔图（蓝色）
for j = 1:length(c_values)
    plot(c_values(j) * ones(1, iterations), bifurcation_data_2(j, :), '.', 'Color', 'b', 'MarkerSize', 1); % 蓝色分岔点
end

xlabel('c');
ylabel('x_n');
title('分岔图');
grid on;
